/***********************************************************************
Compile       : gcc -O3 -o check_solution check_solution.c
Usage         : ./check_solution mode1 mode2 bdflag N L inputfile
                (./check_solution 2 3 1 11 12 Solution.txt)
mode1         : 0 normal ; 1 cycled ; 2 triggled cycled
mode2         : 1 Restriction I ; 2 Restriction II; 3 Restriction I and II
bdflag        : 1 breakdown the designability
N             : number of the nodes
L             : number of the steps

This program is reading the results of the program "design" and checking
whether the solution qualify the additional requirement or not.
./design2 2 0 0 0 0 0 1 11 12 Budding.txt > Solution.txt
./check_solution 2 3 1 11 12 Solution.txt
***********************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

int check_solution (int i, int mode, int N, int state1[N], int state2[N], int state3, int edge[N])
{
	int flag;

	flag = 0;

	if (mode == 1 && state1[i] != state2[i])
	{
		flag = 1;
	}
	else if (mode == 2 && state2[i] == state3)
	{
		flag = 1;
	}
	else if (mode == 3 && (state1[i] != state2[i] || state2[i] == state3))
	{
		flag = 1;
	}

	if (flag == 0)
		return 0;

	int j;

	if (state2[i] == 1)
	{
		/* check inhibition */
		for (j = 0; j < N; j++)
		{
			if (i == j)
				continue;
			if ((state1[j] == 1 || state2[j] == 1) && edge[j] == -1)
				return 1;
		}

		/* check self-activation */
		if (edge[i] == 1)
			return 0;

		/* check activation */
		for (j = 0; j < N; j++)
		{
			if (i == j)
				continue;
			if (state1[j] == 1 && state2[j] == 1 && edge[j] == 1)
				return 0;
		}

		/* check self-inhibition */
		if (edge[i] == -1)
			return 1;

		if (state1[i] == 1)
			return 0;
		else
			return 1;
	}
	else
	{
		/* check inhibition */
		for (j = 0; j < N; j++)
		{
			if (i == j)
				continue;
			if (state1[j] == 1 && state2[j] == 1 && edge[j] == -1)
				return 0;
		}

		/* check self-activation */
		if (edge[i] == 1)
			return 1;

		/* check activation */
		for (j = 0; j < N; j++)
		{
			if (i == j)
				continue;
			if ((state1[j] == 1 || state2[j] == 1) && edge[j] == 1)
				return 1;
		}

		/* check self-inhibition */
		if (edge[i] == -1)
			return 0;

		if (state1[i] == 1)
			return 1;
		else
			return 0;
	}

	return 1;
}

void combine_dist (int N, int L, int count[L], double dist1[N * L], double dist2[N * L], int i)
{
	int j, k;

	if (i == 0)
	{
		for (j = 0; j < L; j++)
		{
			dist2[j] = count[j];
		}
	}
	else
	{
		for (j = 0; j < (i + 1) * L; j++)
		{
			dist2[j] = 0.0;
		}
		for (j = 0; j < L; j++)
		{
			for (k = 0; k < i * L; k++)
			{
				dist2[j + k] += count[j] * dist1[k];
			}
		}
	}
}

int main (int argc, char *argv[])
{
	int mode1, mode2, bdflag, N, L;

	mode1 = atoi (argv[1]);
	mode2 = atoi (argv[2]);
	bdflag = atoi (argv[3]);
	N = atoi (argv[4]);
	L = atoi (argv[5]);

	int pathway[L][N];

	FILE *fp;
	char buf[256];

	int i, j, k;
	double NS, NS2;
	struct solution {
		int NS;
		int *edge;
		int *disqualify;
		int *count;
	} solution[N];

	fp = fopen (argv[6], "r");
	while (1)
	{
		/* read the data of the pathway and the solutions */
		for (i = 0; i < N; i++)
		{
			if (i == 0 && mode1 == 2)
				continue;
			if (fgets (buf, 256, fp) == NULL)
				break;
			sscanf (buf + 39, "%d", &solution[i].NS);
			solution[i].edge = calloc (solution[i].NS, sizeof (int) * N);
			solution[i].disqualify = calloc (solution[i].NS, sizeof (int));
			solution[i].count = calloc (L - 1, sizeof (int));
			for (j = 0; j < solution[i].NS; j++)
			{
				fgets (buf, 256, fp);
				for (k = 0; k < N; k++)
				{
					if (buf[10 + 2 * k] == 'n')
						solution[i].edge[j * N + k] = 0;
					else if (buf[10 + 2 * k] == 'g')
						solution[i].edge[j * N + k] = 1;
					else if (buf[10 + 2 * k] == 'r')
						solution[i].edge[j * N + k] = -1;
				}
			}
		}
		if (i != N) break;
        
		fgets (buf, 256, fp);
        
		for (i = 0; i < L; i++)
		{
			for (j = 0; j < N; j++)
			{
				fscanf (fp, "%d", &pathway[i][j]);
			}
			fgets (buf, 256, fp);
		}
		/* end of read */
        
		/* check the solutions */
		for (i = 0; i < N; i++)
		{
			if (i == 0 && mode1 == 2)
				continue;
        
			for (j = 0; j < solution[i].NS; j++)
			{
				if (mode1 != 2)
				{
					solution[i].disqualify[j] = check_solution (i, mode2, N, pathway[0], pathway[1], pathway[2][i], solution[i].edge + j * N);
				}
				for (k = 1; k < L - 2; k++)
				{
					solution[i].disqualify[j] += check_solution (i, mode2, N, pathway[k], pathway[k + 1], pathway[k + 2][i], solution[i].edge + j * N);
				}
				if (mode1 == 0)
				{
					solution[i].disqualify[j] += check_solution (i, mode2, N, pathway[k], pathway[k + 1], pathway[k + 1][i], solution[i].edge + j * N);
				}
				else
				{
					solution[i].disqualify[j] += check_solution (i, mode2, N, pathway[k], pathway[k + 1], pathway[1][i], solution[i].edge + j * N);
				}
        
				/*
				printf ("Node %2d :", i);
				for (k = 0; k < N; k++)
				{
					if (solution[i].edge[j * N + k] == 0)
						printf (" n");
					else if (solution[i].edge[j * N + k] == 1)
						printf (" g");
					else if (solution[i].edge[j * N + k] == -1)
						printf (" r");
				}
				printf (" %2d\n", solution[i].disqualify[j]);
				*/
			}
        
			for (j = 0; j < L - 1; j++)
			{
				solution[i].count[j] = 0;
			}
			for (j = 0; j < solution[i].NS; j++)
			{
				solution[i].count[solution[i].disqualify[j]] ++;
			}
		}
		/* end of check */
        
		/* calculate the distribution */
		double dist1[N * (L - 1)], dist2[N * (L - 1)];
		for (i = 0; i < N; i++)
		{
			if (i == 0 && mode1 == 2)
				continue;
        
			if (mode1 == 2)
				combine_dist (N, L - 1, solution[i].count, dist1, dist2, i - 1);
			else
				combine_dist (N, L - 1, solution[i].count, dist1, dist2, i);
        
			memcpy (dist1, dist2, sizeof (double) * N * (L - 1));
		}
        
		if (bdflag == 1)
		{
			for (i = 0; i < N * (L - 1); i++)
			{
				if (dist2[i] > 0.9)
				{
					printf ("Designability for %2d : %lg\n", i, dist2[i]);
				}
			}
		}
		/* end of calculate */

		NS = 1.0;
		for (i = 0; i < N; i++)
		{
			if (i == 0 && mode1 == 2)
				continue;

			NS *= solution[i].NS;
		}

		NS2 = 0.0;
		for (i = 0; i < N * (L - 1); i++)
		{
			NS2 += dist2[i];
			if (NS2 * 2.0 > NS)
				break;
		}
		printf ("Total Designability : %lg\tMedian : %2d\n", NS, i);
	}

	fclose (fp);
}
